Molecular dynamics with quantum heat baths: results for nanoribbons and nanotubes 
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A generalized Langevin equation with quantum baths (QMD) for thermal transport is derived 
with the help of nonequilibrium Green's function (NEGF) formulation. The exact relationship of 
the quasi-classical approximation to NEGF is demonstrated using Feynman diagrams of the nonlin- 
ear self energies. To leading order, the retarded self energies agree, but QMD and NEGF differ in 
lesser/greater self energies. An implementation for general systems using Cholesky decomposition of 
the correlated noises is discussed. Some means of stabilizing the dynamics are given. Thermal con- 
ductance results for graphene strips under strain and temperature dependence of carbon nanotubes 
are presented. The "quantum correction" method is critically examined. 
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I. INTRODUCTION 

Molecular dynamics 1 (MD) has been used as one of 
the most important simulation tools to study a variety 
of problems from structure to dynamics. In particular, 
MD is routinely used for thermal transport problems^. 
MD is versatile and can handle with ease any form of 
classical interaction forces. The computer implementa- 
tion of the algorithms is straightforward in most cases. 
However, there is one essential drawback in MD - it is 
purely classical, thus it is unable to predict quantum be- 
havior. Of course, in situations like high temperatures 
and atoms other than hydrogen, MD gives good approx- 
imation. This is not the case for very small systems like 
nanostructures at low temperatures. For lattice vibra- 
tions, the relevant temperature scale is the Debye tem- 
perature, which is quite high for carbon-based materials. 
Thus, even the room temperature of 300 K is already 
considered a low temperature. One uses MD anyway in 
cases where quantum effects might be important, for lack 
of better alternative approaches. 

Recently, it was proposed that MD can be augmented 
with a quantum heat bath to at least partially take into 
account the quantum effect*. We'll refer to this new 
molecular dynamics as QMD. The proposed generalized 
Langevin dynamics with correlated noise obtained ac- 
cording to Bosc distribution has the important features 
that it gives correct results in two special limits, the low- 
temperature ballistic limit and high-temperature diffu- 
sive limit. It is one of the very few methods that bal- 
listic to diffusive transport can be studied in a single 
unified framework. The QMD should be most accurate 
for systems with strong center- lead couplings. The non- 
Markovian heat baths have an additional advantage in 
comparison with the usual Langevin or Nose-Hoover heat 
baths in that the baths (the leads) and the systems can 
be connected seamlessly without thermal boundary re- 
sistance. 

In this paper, we follow up the work of ref. |4| to give 
further details and calculations on large systems. Wc 
first discuss the implementation of the colored noises with 



several degrees of freedom. We then discuss the problem 
of instabilities we encountered in carrying out the sim- 
ulation. Various ways of overcoming this difficulty are 
suggested and tested. We analyze the proposed dynam- 
ics and compare it with the exact nonequilibrium Green's 
function (NEGF) method^ in terms of Feynman diagrams 
and nonlinear self-energies. Here we show that in the bal- 
listic case, the generalized Langevin dynamics and NEGF 
are completely equivalent. Using this analysis, it is also 
clear that at high temperatures NEGF and Langevin dy- 
namics should agree for nonlinear systems. We then re- 
port some of the simulation results on nanoribbons (with 
periodic boundary condition in the transverse direction) 
and nanotubes. We comment on one popular method 
of "quantum correction" and point out its shortcomings 
and inconsistencies. We conclude in the last section. 



II. GENERALIZED LANGEVIN DYNAMICS 
AND IMPLEMENTATION DETAILS 

The derivation of the generalized Langevin equation 
for junction systems was given in ref. 0, see also refs. 
and |7|. Here we present a much faster derivation using 
the results of NEGF, following the notations introduced 
in ref.0. The starting point is the set of quantum Heisen- 
berg equations of motion for the leads and center, 
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u a is a vector of displacements in region a away from 
the equilibrium positions, multiplied by the square root 
of mass of the atoms. The leads and the coupling be- 
tween the leads and center are linear; while the force in 
the center, F c (u c ) = -K c u c + F n , is arbitrary. We 
eliminate the lead variables by solving the second equa- 
tion and substituting it back into the first equation. The 
general solution for the left lead is 
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where g r L {t) is the retarded Green's function of a free 
left lead with the spring constant K L , satisfying g r L (t) + 
K L g r L {t) = -8(t)I, g r L (t) = for t < 0. Its Fourier 
transform is given by [(w + ir/) 2 — K L ]~ 1 , where 77 — > + 
is an infinitesimal positive quantity, u^it) satisfies the 
homogeneous equation of the free left lead: 



K u. 
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g r L (t) and Uq (t) are associated with the "free" lead in the 
sense that leads and center are decoupled, as if V CL = 0. 
This is consistent with an adiabatic switch-on of the lead- 
center couplings. The right lead equations are similar. 
Substituting Eq.([3]) into Eq. ([!]), we obtain 
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where E r = 17 L + Y7 R , TT a = V Ca g r a (t)V aC ', is the self- 
energy of the leads, and the noise is defined by = 
-V Ca u%(t), £ = 

The most important characterization of the system is 
the properties of the noises. This is fixed by assuming 
that the leads are in respective thermal equilibrium at 
temperature Tl and Tr. It is obvious that for a set of 
coupled harmonic oscillators, there is no thermal expan- 
sion effect, (iig(i)} = 0, thus (£ a ) = 0. The correlation 
function of the noise is 
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where the superscript T stands for matrix transpose. We 
have used the definition of greater Green's function and 
self-energy of the free left lead 5 . We assume that the 
noises of the left lead and right lead are independent. 
Since the noises £ a {t) are quantum operators, they do 
not commute in general. In fact, the correlation in the 
reverse order is given by the lesser self-energy: 
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Equation ([5]) together with the noise correlations 
Eq. © and (0 is equivalent to NEGF approach. For 
the quantum Langevin equation, it is not sufficient to 
completely characterize the solution by just the first and 
second moments of the noises. We need the complete 
set of n-point correlators (£(ti)£(t2) ■ ■ • £(*n))i which is 
in principle calculable from the equilibrium properties of 
the lead subsystem 8 . It is very difficult to solve the dy- 
namics unless the nonlinear force F n is zero. Thus, for 
computer simulation in the quantum molecular dynam- 
ics approach, we have replaced all operators by numbers 
and a symmetrized noise, ihh[^ < (t) + S > (t)] = £7i£(t), 
is used. This is known as quasi-classical approximation 
in the literature 9,10 . 



A. Implementation 



The formula for the noise spectrum of the left lead is 



F[oj] = iKEi[w] = h 
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where f(u) = \/[ e hu ^ k B T L) _ 1] { s the Bose distribution 
function, and T L [u>] = i (Y, r L [lu] — M) ■ The right lead is 
analogous. The surface Green's functions g r a are obtained 
using an iterative methodi^iii To generate the multivari- 
ate gaussian distribution with an arbitrary correlation 
matrix, we use the algorithm discussed in ref. [l2l . That 
is, we do Z — cX, where X is a complex vector following 
standard uncorrelated gaussian with unit variance, one 
for each discretized frequency, while cc T — F[u>], and c is 
a lower triangular real matrix, c is obtained by Cholesky 
factorization?^ from a Lapack routine dpotrf ( ) . The 
Cholesky decomposition is performed only once. The 
frequency array of lower triangular matrices c is stored. 
The Fourier transform of Z gives the noise in time do- 
main, which is obtained using a fast Fourier transform 
algorithm. Further details are given in ref. EL 



B. Overcoming instability 

We have implemented the second generation reactive 
empirical bond order (REBO) Brenner potential^ for 
carbon with the special restriction that the coordination 
numbers are always three. This is valid for carbon nan- 
otubes and graphene sheets with small vibrations in ther- 
mal transport. We found that a naive implementation of 
the QMD in higher dimensions unstable. The atoms close 
to the leads have a tendency to run away from the poten- 
tial minima and go to infinity. Several ways were tried 
to stabilize the system. 

(1) Instead of integrating over the coordinates u (t) 
in the memory kernel, we can perform an integration by 
part, and consider integrating over velocity. This form 
of the generalized Langevin equation resembles more of 
the standard Langevin equation of velocity damping, but 
there will be an extra force constant term, as follows: 

u c = F c + \r{0)u c -J Y{t-t')u c (t')dt' + f, (9) 

where T(t) is defined by Eq. (fl3|) below. We introduce a 
parameter A, which should take the value 1, but using a 
smaller value can stabilize the system. However, A / 1 
introduces boundary resistance. 

(2) We scale up the force constants of the leads by a 
factor of /. This broadens the lead spectra to be closer 
to white noise, thus better damping. 

(3) We add an additional onsite force on each atom, 
with a linear force constant if onsite, as well as a small 
u nonlinear force. This breaks the translational invari- 
ance so that the atoms are fixed near their equilibrium 
positions. 
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FIG. 1: The surface density of states D s (u>) vs. frequency 
for a (n, 2) zigzag graphene strip with (1,2) of 8 atoms as a 
repeating unit cell. The delta peaks are located at 566, 734, 
1208, 1259, 1287, and 1632 cm" 1 . The rest of the peaks do 
not diverge as r\ — > 0. 
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FIG. 2: (color online). Normalized vibration amplitudes 
vs. reduced coordinate x of each carbon atom. From (a) to (f) 
are six edge modes in (10, 2) graphene strip (blue solid line) 
and (20, 2) graphene strip (red dotted line). The frequency 
u for each mode given in the figure is in cm" 1 . 



(4) We smooth the noise spectrum by choosing a small 
number of points, say 100 sampling points in frequency. 
We add an artificial damping, e~ et to Eq. (fT3|l . 

(5) We implemented three algorithms: velocity Ver- 
let, fourth order Runge-Kutta, and an implicit two-stage 
fourth-order Runge-Kutta^. 

Not all of the measures are effective. We feel perhaps 
the most important point is (4). The extra parameters 
A, /, -Konsite, and e will be stated when discussing the 
results. 



III. DELTA-PEAK SINGULARITIES IN LEAD 
SELF-ENERGY 

For a one-dimensional (ID) harmonic chain with a uni- 
form spring constant K, the lead self-energy is given by 
E r [w] = -if A, where A satisfies K A + (u> + ir,) 2 -2K + 
K/\ = and |A| < 1. Both the real part and imagi- 
nary part are smooth functions of the angular frequency 
ui. However, this is not true in general. For sufficiently 
complex leads, we find (5-function-like peaks on an other- 
wise smooth background, see Fig. [TJ What is plotted in 
Fig. [1] is the surface density of states defined according 
to 

D s (ui) = -— ImTr<?£[u;]. (10) 

7T 

The above formula gives the bulk phonon density of 
states if g r is replaced by the central part Green's func- 
tion G r . The peaks in Fig. Q] are not numerical arti- 
facts, but real singularities in the semi-infinite lead sur- 
face Green's function or self-energy. If we omit them, the 



identity (a special case of the Kramers-Kronig relation) 

J-oo 7T W 

will be violated. These sparks are indeed 5 functions. As 
the small quantity r\ in (lu + irj) 2 decreases, the peaks 
become higher and narrower, but the integral in a fixed 
interval around the peak remains constant. These peaks 
are not related to the van Hove singularities of bulk sys- 
tem density of states, as the locations of the peaks are 
not at these associated with zero group velocities. The 
singularities of the self-energy can be approximated as 
proportional to 

-wP — iir5(uj~uJo), (12) 

lu — luq + irj lu — luq 

where P stands for principle value. 

To interpret these peaks, we did a calculation of the 
vibrational eigenmodes for a finite system. We find 
localized modes with frequencies matching that of the 
delta peaks in D s (lu). We use "General Utility Lattice 
Program" (GULP)^ to calculate the phonon modes of 
graphene strip, with fixed boundary condition in x di- 
rection and periodic boundary condition in y direction. 
These boundary conditions are the same as that in the 
calculation of lead surface Green's function. Six local- 
ized modes are found in (n, 2) graphene strips. Fig. [2] 
shows the normalized vibrational amplitude of each atom 
in all six localized modes of (10,2) (blue solid) and (20, 
2) graphene strip (red dotted). In these modes the vi- 
brational amplitude decreases exponentially to zero from 
edge into the center. The frequency of these local- 
ized modes are 561.5, 741.3, 1193.5, 1264.5, 1279.6, and 
1636.0 cm" 1 . These values match very well to that of the 
delta peaks in D s (uj) shown in Fig. [TJ These frequency 
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values are the same in (10, 2) and (20, 2) graphene strip. 
After a critical distance L c from the edge, the vibration 
amplitude decreases to zero. L c is the same in (10, 2) 
and (20, 2) graphene strip. So these localized modes 
are relatively more 'localized' in longer graphene strip 
as shown in Fig. O Because of their localizing prop- 
erty, these modes are important in thermal transport. 
There are localized modes both at the edges of leads and 
edges of center. They have opposite effect on the ther- 
mal conductance. (1) Localized modes at the edges of 
leads are beneficial for thermal conductance. Because 
in these modes, atoms at the edges of leads have very 
large vibration amplitude. As a result, thermal energy 
can transport from leads into center more easily. (2) 
Localized modes at the edges of the center have negative 
effect on thermal conductance. In these modes, only out- 
side atoms have large vibrational amplitudes while inside 
atoms have small vibration amplitudes or even do not vi- 
brate at all. So thermal energy are also localized at the 
edges, making it difficult to be transported from one end 
to the other end. The effects of (1) and (2) are such that 
the net result is equivalent to a perfect periodic system 
without boundary resistance. 



The localized modes are a consequence of dividing the 
infinite system abruptly and artificially into leads and 
center. Implementing these delta-peak singularities in 
a QMD simulation is impossible, since these modes are 
not decaying in time for the real-time self-energy £ r (i). 
Thus, we are forced to remove these peaks from the imag- 
inary part of E r , and reconstruct a real part using the 
Hilbert transform from the imaginary part with the delta 
peaks removed. The damping kernel for the QMD dy- 
namics is computed from (for t > 0) 



IV. DIFFERENCE BETWEEN QMD AND 
NEGF: A FEYNMAN DIAGRAMMATIC 
ANALYSIS 

In this section, we give an analysis of the difference 
between fully quantum-mechanical NEGF and the quasi- 
classical generalized Langevin dynamics. A similar result 
was presented in ref. [171 briefly for the case of electron- 
phonon interaction. The starting point is a formal so- 
lution of Eq. ([5]) with the quasi-classical approximation 
and a symmetrized correlation matrix for the noise: 



u(t) 



G r (t : t')[^t')+F n (t')]dt\ (14) 



where we have omitted the superscript C on u for sim- 
plicity, G r is the retarded Green's function of the central 
region for the ballistic system (when F n = 0). We have 
also left out a possible term satisfying a homogeneous 
equation [Eq. ([5]) when £ = F n = 0] and depending on 
the initial conditions. Physically, such term should be 
damped out. Provided that the central part is finite, such 
term should not be there and this is consistent with the 
fact that the final results are independent of the initial 
distribution of the central part in steady states. 

We consider the expansion of the nonlinear force of the 
form 



{F n )i 



^ TijkiUjUkUi, 
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(15) 



where T^-fe and Tijki are completely symmetric with re- 
spect to the permutation of the indices. From repeated 
substitution of Eq. (fT4"|) back into itself, we can see that 
u(t) is expressed as polynomials of G r and £. The corre- 
lation functions of u can then be calculated using the fact 
that the noise is gaussian, and Wick's theorem applies. 
It is advantageous to define two types of (quasi-classical) 
Green's functions, as 



r(t) 



Ys r {t')dt' 



duj ImE [lu] 
c 

7r — lj + ie 



(13) 

In practice, the removal of the peaks is done by choosing 
a small r\ (ss 10~ 8 ). Since the sampling of uj is at a 
finite spacing, typically with about 10 2 points, we almost 
always miss the peaks if rj is small. 



In calculating the ballistic transmission through Caroli 
formula, the omission of the delta peaks at a set of points 
of measure zero has no consequence. However, the exis- 
tence of the singularities is also reflected through the real 
part of the self-energy. If the real part uses the Hilbert 
transformed version with the delta peaks omitted, the 
transmission coefficient T[oj] will not be flat steps as ex- 
pected for a perfect periodic system. Thus, removing the 
delta peaks consistently means we are using a lead that 
is modified from the original one. 



-( u {t)u T {t')) = g n (t,t') 



(16) 

g r n (t,t")x L (t", t')dt", (l?) 



The energy current is calculated by the amount of de- 
crease of energy in the left (or right) lead: 



(18) 



We can replace u L with the solution, Eq. . Going into 
the Fourier space and some algebraic manipulation, we 
can write 



II = - 



— hu Tr(<?; MS L [w]+a n [w]EJ M) , (19) 



where the Fourier transform is defined in the usual 
way, e.g., Q[oj] = f_°° Q(t)e luJt dt. The above equa- 
tion has the same form as the NEGF one, provided 
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FIG. 3: The top two lines are for the quasi-classical self- 
energies and £ n ; the next two lines are the corresponding 
NEGF results. A line without arrow represents G. A line 
with an arrow represents G r when read following the arrow, 
or G a when read against the direction of arrow. A line with 
single sided arrow represents G > when following the arrow 
and G < when read against the sense of arrow. 



that we can identify the quasi-classical Green's functions 
defined in Eq. (|16|) and (JTTJ) with the quantum ones, 
G = i(G> + G<) and G r . [It looks slightly different 
from the expression of Eq. (5) in ref. 18, where only G < 
appears. There is an error in that paper. One should 
take the real part of that expression, or add its complex 
conjugate. By doing this, we obtain a symmetrized ex- 
pression with respect to G < and G > .] 

We compare Q r n and Q n with their fully quantum- 
mechanical counterpart G r n and G n through the nonlin- 
ear self-energies. It is not known if a Dyson equation 
for Q r n is still valid in the sense that the self-energy con- 
tains only irreducible graphs, but we simply 'define' the 
retarded nonlinear self-energy through 



EI = G r 



and similarly define E„ by 



g n = Ql (£ + t n ) G a n . 



(20) 



(21) 



To simplify the representation of the diagrams, we have 
used the fact that G r jt {t,t') = Gfj(t',t), Gf t (t,t') = 
Gy(t',t), and the Keldysh relation in frequency domain 
G = G 7 "SG a . In Fig. [3] we give the lowest-order dia- 
grams of the two types of self-energies and contrast with 
the NEGF results. Numerous cross terms involving prod- 
ucts of Tijk and Tiju are not shown. The NEGF results 
are obtained in ref. Qj^ (Fig. 3) for the contour ordered 
version, here we have separated out explicitly for Y7 n and 
E n = ±(E> + E<). It is clear from Eq. and {21]) 
that, when the nonlinear couplings T^k and Tijki are 
zero, we have Q r n = G r and Q n = G. Thus, for ballis- 
tic systems, NEGF and quasi-classical MD agree exactly. 
To leading order in the non-linear couplings [0(T^ k ) and 




FIG. 4: (color online), 
graphene strip with (n, m) 
the translational period. 



The structure for an armchair 
= (4,2). The box in red line is 



0(Tijki)] the retarded nonlinear self-energies agree. The 
difference starts only at a higher order. The self-energies 
E„ disagree even at the lowest order. The NEGF and 
quasi-classical diagrams become the same if we take the 
"classical limit" Ti — + with a new definition of classical 
Green's functions G r -> G r d and TiG> w hG< -> G d . 
In this limit, the distinction between G > and G < disap- 
pears. The extra diagrams go to zero because they are 
high orders in h. 



V. TESTING RUNS AND COMPARISON WITH 
NEGF 

Fig. [1] is the configuration of a system in the simula- 
tion. There are four atoms in the translational period. 
A pair of numbers (n, m) are introduced to denote the 
number of periods in the horizontal and vertical direc- 
tions. They should not be confused with the chirality 
indices of the nanotubes. This figure shows the particu- 
lar case of armchair graphene strip with (n, m) = (4, 2). 
For zigzag configurations, the unit cell is rotated by 90 
degrees. In the vertical direction, a periodical boundary 
condition is applied. In the simulation box, atoms in the 
left-most columns labelled 0-7 are fixed left lead, atoms 
in the right-most columns labelled 28-35 are the fixed 
right lead, and the heat baths are applied to the columns 
close to them. The temperature of the leads are set ac- 
cording to T L = T(l + a), and T R = T(l - a). The ther- 
mal conductance is computed from a = Il/(Tl — Tr). 

Test runs, shown in Fig. [5l with parameters a = 0.4, 
A = 0.6, / = 1.2, K onsitc = 0.01 (eV/(A 2 u)), and 
e = 0.001 (10 14 Hz) with the geometry of Fig.H demon- 
strate that QMD implemented by the velocity Verlet and 
fourth order Runge-Kutta gives the correct results in 
comparison with ballistic NEGF. For a system of such 
small sizes, the conductance behaves ballistically. The 
one implemented by the velocity Verlet agrees very well 
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FIG. 5: (color online). A comparison of temperature depen- 
dence of the thermal conductance for an armchair graphene 
strip with (n,m) = (4,2), solid line: NEGF, circle: QMD 
with velocity Verlet, square: QMD with fourth order Runge- 
Kutta. (a) The A dependence for the same system at 300 K 
for QMD with velocity Verlet (circle) and NEGF (solid line). 



with the NEGF result in the low-temperature regime. 
Other implementation methods, like an implicit two- 
stage fourth order Runge-Kutta, also turn out to give 
similar results. Thus, the results are rather insensitive to 
the integration algorithms used. This suggests the suc- 
cess of simulating quantum transport not only for the 
one-dimensional quartic onsite model^, but also for the 
large systems. Due to the artificial parameters added in 
order to overcome the instability, the thermal conduc- 
tance obtained was slightly higher than the ballistic one 
in high-temperature regime. We note that the param- 
eters A, /, and K ons [ te are incorporated in the NEGF 
calculation, the effect of e is not taken care correctly in 
NEGF. This may explain the discrepancy at high tem- 
peratures. We further analyze the A dependence of the 
thermal conductance for the (4, 2) graphene strip: the 
inset Fig. 02a) represents the room temperature (300 K) 
results, where the thermal conductance exhibits linear 
dependence on A. The conductance reduces by about half 
when A is reduced from I to 0. Besides A, other param- 
eters also have their own impacts, for instance, smaller 
e lowers the effect of the artificial damping, but requires 
much larger integration domain and therefore brings the 
risk of truncating the spectrum and providing the wrong 
self energy. The conductance is independent of e if it is 
in the range 0.001 to 0.02. 
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FIG. 6: (color online). The dependence of the thermal con- 
ductance on the length of the system at 300 K in double loga- 
rithmic scale. Phonon transport changes gradually from bal- 
listic to diffusive with increasing length of the system. 



0.22 




FIG. 7: (color online). The effect of strain on the thermal 
conductance of the graphene. The dotted line (in red) is guide 
to the eye. 



VI. RESULTS OF NANORIBBON UNDER 
TENSION 

In the simulation, typical MD steps are 10 5 of 0.5 fs, 
which is long enough to obtain converged thermal con- 
ductance. The stabilizing factor is A = 0.6 with an onsite 
force constant for all the atoms K ons i tc = 0.01 eV/(A 2 u), 
and other parameters / = 1.2 and e — 10 -3 (10 14 Hz). 

Fig. [6] shows the dependence of the thermal conduc- 
tance on the length (L) of the system for (n, 2) zigzag 
graphene strip at 300 K. There are three regions labeled 
(I), (II) and (III) in the figure respectively for L in [10, 
100], [100, 600] and [600, 1400] A. In the very short 
length region (I), thermal transport should be in the bal- 
listic regime, where thermal conductance is a constant 
independent of the length of the system. But here the 
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FIG. 8: (color online). Phonon density of state (DOS) for 
graphene under elongation strain (a), and compression strain 
(b). 



thermal conductance exhibits decreasing behavior. Ac- 
tually, this decreasing behavior is mainly attributed to 
the localized modes at the edge of the center region. As 
shown in Fig. [21 these modes will be more localized in 
longer graphene strips. So they make little contribution 
to the thermal transport in short graphene strips due to 
their localizing property, and even this little contribution 
is further reduced with increasing length. That is the 
reason for the decrease of thermal conductance in this 
length region. After L w 100 A, these localized modes 
are fully localized, so they do not contribute to the ther- 
mal conductance anymore. In region (II), the thermal 
transport is in the ballistic regime, where the thermal 
conductance is more or less a constant. In region (III), 
this curve decreases as L increases, indicating cross-over 
to the diffusive thermal transport. This length scale is 
consistent with previous theoretical results.— Thermal 
conductance in this region can be fitted by a power func- 
tion a — 6.1 £~ - 66 (dotted line). So the thermal con- 
ductivity is proportional to the length as L 34 . This 
exponent 0.34 agrees with previous results on nanotubes 
and other quasi-lD systems . 21 i 22 i 23 

The associated values of thermal conductivity, k = 
aL/S, where L is length and S is cross-section area, is 
too small. The smallness is attributed to the bound- 
ary resistance caused mainly by A ^ 1 and the omis- 
sion of the delta-peak lead self energies. For a perfect 
ID system of (oo, 2) the conductance with the Brenner 
potential is 0.72 nW/K from a ballistic NEGF calcula- 
tion. If the leads are replaced by those of omitting the 
delta-peaks as discussed in Sec. Ill, the NEGF (4,2) sys- 
tem result that is consistent with our simulation setup 
is reduced to 0.19 nW/K. This is quite close to, but still 
some discrepancy with, QMD result. These may be due 
to nonlincarity and other unexplained systematic errors. 

In Fig. [71 the effect of the strain on the thermal con- 
ductance of graphene is displayed for a system of zigzag 



(4, 2). To mimic the experimental condition ; 24 ' 2 ^ 26 the 
strain is introduced in two steps. Firstly, the strain is 
generated to the whole graphene system in Fig. [4j Sec- 
ondly, atoms in the center are fully relaxed with left and 
right leads fixed. And we then do the MD simulation on 
this optimized graphene system. We find that the ther- 
mal conductance increases with increasing elongation on 
graphene. But compression on graphene does not change 
the value of the thermal conductance appreciably. 

To understand this strain effect on the thermal conduc- 
tance, we study the density of state (DOS) of the phonons 
in Fig. [H The DOS is calculated from the Brenner em- 
pirical potential as implemented in GULP for (4, 2) ge- 
ometry with fixed boundary conditions in x direction, 
and periodically extended in y direction. We use GULP 
to do optimization for the strained graphene with two 
leads fixed firstly, and then calculate the DOS of this re- 
laxed system. As shown in Fig. [5] (a), the high frequency 
Raman active mode (G mode) around 1600 cm -1 shows 
obvious red-shift under extension strain, which agrees 
with the recent experimental results ; 24 ! 25 Furthermore, 
Fig. [8] (b) predicts the blue-shift of the G mode with 
compression strain. 

For thermal conductance of the graphene at room 
temperature, the phonon modes with frequency about 
200 cm -1 are important. We can see two significant 
modes (1 and 2) in this frequency region in Fig. [5] When 
the graphene is elongated, both modes 1 and 2 are blue- 
shifted [Fig. [8] (a)], which results in increasing of thermal 
conductance. However, if the graphene is compressed, 
modes 1 and 2 shift in opposite directions [shown in 
Fig. [8] (b)]. As a result, the contribution of these two 
modes to the thermal conductance cancels with each 
other. That is the reason for the almost unchanged value 
of the thermal conductance under compression. 



VII. RESULTS ON NANOTUBES 

Fig. [9] shows our simulation results on zigzag carbon 
nanotubes of chirality (5,0) with different lengths by us- 
ing the same parameters as that of the previous sec- 
tions. Each data point typically takes about 48 hours 
on an AMD Opteron CPU. The thermal conductivity 
is computed according to k = aL/S, where L is the 
length of the sample, and assuming a cross-section area of 
S = 12 A 2 . Both the thermal conductance and thermal 
conductivity monotonically increase with the tempera- 
ture in the low-temperature regime, which agrees with 
the available experimental data and demonstrates the 
ability of QMD to predict the quantum effect in this 
regime. This is completely neglected in the classical 
MD approaches . 22 ' 27 i 28 ' 29 i 30 For nanotubes with 12.8 nm 
(30,5) and 25.6 nm (60,5) length, as the temperature 
increases, the thermal conductance and corresponding 
thermal conductivity start to drop at 850K, this decre- 
ment is consistent with the classical prediction, which 
indicates that the quantum correction becomes much 
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FIG. 9: (color online). Temperature dependence of the ther- 
mal conductance and thermal conductivity for zigzag carbon 
nanotubes with (n,m) — (10,5) (triangle), (30,5) (square), 
and (60, 5) (circle). 



smaller. Yet, such decline has not been observed in the 
shorter case with a length of 4.26 nm. This difference 
shows that a transition from ballistic to diffusive happens 
when the length gets longer. The thermal conductance 
decreases but conductivity still increases with nanotubes 
length; we are still in transition to diffusive regime^ 
However, the values at high temperatures are compara- 
ble to previous MD results. 



VIII. A CRITIQUE TO THE QUANTUM 
CORRECTION METHOD 



The quantum correction method was first sug- 
gested in refs. I3lll32lf33l and used by a number of 
researcher s 30 ' 34 ! 35 ' 36 without carefully examining its va- 
lidity. From the simple kinetic theory of thermal trans- 
port coefficient, the thermal conductivity can be written 
as k = l/(3V)J2k c k v klk where Ck is the heat capacity 
of a mode fc, and Vk and Ik are the associated phonon 
group velocity and mean free path, V is the volume of 
the system. Provided that the phonon velocity Vk and 
mean free path Ik of mode k are approximately indepen- 
dent of k, we can argue that the quantum conductivity 
is scaled down by the quantum heat capacity from the 
classical value. In quantum correction, the temperature 
is also redefined such that the classical kinetic energy is 
equated with the corresponding quantum kinetic energy 
of a harmonic lattice. Here it is not clear whether the 
zero-point motion should be included or not . 30 ' 32 ' 37 

To a large extent, a constant phonon velocity is a 
good approximation. However, it is well-known that the 
phonon mean free path is strongly dependent on the fre- 



quencies, e.g., in Klemens' theory for umklapp process 3 -^, 
I ~ uj~ 2 . To what extent the quantum correction works 
is questionable. There is one special case that we can 
answer this question definitely, although this is for con- 
ductance, not conductivity. Let us consider a ID har- 
monic chain and compute its conductance exactly and 
compare it with a classical dynamics with quantum cor- 
rection. The correct answer for the thermal conductance 
is given by the Landauer formula, 



CQM 



— fruj T \uj] -J- 



(22) 



where the transmission T[uj] is one for a uniform chain, 
wmax = y/4K/m is the maximum frequency for a chain 
with spring constant K and mass m. c(w) = Hcudf/dT is 
the heat capacity of the mode at frequency u>. The cor- 
responding classical value is obtained by approximating 
the Bose distribution function with / « k B T '/ 'Quo). This 
gives the correct classical value of conductance as 



ocl = — - — k B . 

Zir 



(23) 



Now we consider quantum correction to Eq. ([23)1 . The 
total quantum heat capacity of a ID harmonic chain is 



C = J2 



Ck 



L 



c{u>) duj 

v(uj) 7T 



(24) 



where v(lu) — duo/dk — (a/2)y / w^ IAX — to 2 is the phonon 
group velocity. The classical value is Nk B — Lk B /a, a is 
lattice constant. 

According to the quantum correction scheme, the re- 
sult from a classical dynamics is corrected by multiplying 
the classical value by the ratio of quantum to classical 
heat capacity, given 

" MAX a^MAX . ,dw 

TTV(UJ) Zir 



C f 
Nk B Jo 



This does not agree with the correct quantum result of 
Eq. (f2"2"]l . There is no need to shift the classical tem- 
perature as <7cl is independent of the temperature. The 
heat capacity at frequency u> is weighted differently in 
two cases. Even if the group velocity can be approxi- 
mated by a constant by v(0) = a^J K/m, valid at very 
low temperatures, the two results still differ by a factor 
of tt/2. 



IX. CONCLUSION 

We have presented a quick derivation of the general- 
ized Langevin equation, emphasizing its connection with 
NEGF. The inputs to run the Langevin dynamics can 
be calculated in the standard way from a NEGF phonon 
transport calculation. The implementation details are 
given, such as the generation of colored noise vector £. 
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We found quite generically that the lead self-energies 
contain delta-function peaks for quasi-one-dimensional 
systems. These delta peaks represent surface or edge 
modes. This complicates the molecular dynamics sim- 
ulations. These delta peaks in the spectra have to be 
removed in order to obtain a stable simulation. We hope 
that the instability is specific to the systems of quasi-lD 
carbon graphene strips or nanotubes. If the leads are 
modeled as bulk 3D systems, the noise spectra should 
be more smooth, and should produce a stable dynamics. 
The quasi-classical approximation which results to the 
generalized Langevin equation is analyzed using Feyn- 
man diagrams and its results are compared with NEGF. 
It is found that, to lead order, the nonlinear retarded 
self-energy agrees with NEGF, while £„ does not, mainly 
due to the fact that QMD cannot distinguish between G < 
and G > . As a by-product, we see easily that QMD and 
NEGF agree for linear systems. QMD also gives the cor- 



rect classical limit. We presented test runs and compared 
with NEGF for the thermal conductance. Long (n, 2) 
graphene strips are simulated to study the crossover from 
ballistic transport towards diffusive transport. Effect of 
strain is also studied. The results of carbon nanotubes 
are also presented. Our simulations are one of the first 
examples of the QMD on realistic systems. Finally, the 
quantum correction method is critically examined. 
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